Detailed Look
A more detailed look at spatial interpolation, integration through time, and I/O. For additional documentation e.g. see 1, 2, 3, 4. Here we illustrate a few things in more detail:
- reading velocities from file.
- gridded velocity output (Udata, Vdata)
- pre-computed trajectory output (
float_traj*data)
- interpolating
U,Vfrom gridded output to individual locations- compared with
u,vfromfloat_traj*data
- compared with
- computing trajectories (location v time) using
OrdinaryDiffEq.jl- compared with
x(t),y(t)fromfloat_traj*data
- compared with
1. Import Software
using IndividualDisplacements, OrdinaryDiffEq, DataFrames, MITgcmTools
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_plots.jl"))
include(joinpath(p,"../examples/example123.jl"))
include(joinpath(p,"../examples/helper_functions.jl"));2. Read Trajectory Output
from MITgcm/pkg/flt
dirIn=IndividualDisplacements.flt_example
prec=Float32
df=read_flt(dirIn,prec);
plt=plot_paths(df,300,100000.0)
3. Read Gridded Variables
using MeshArrays.jl and e.g. a NamedTyple
𝑃,Γ=example2_setup();4. Visualize Velocity Fields
plt=heatmap(Γ["mskW"][1,1].*𝑃.u0,title="U at the start")
plt=heatmap(Γ["mskW"][1,1].*𝑃.u1-𝑃.u0,title="U end - U start")5. Visualize Trajectories
(select one trajectory)
tmp=df[df.ID .== 200, :]
tmp[1:4,:]| ID | time | lon | lat | dep | i | j | k | etaN | uVel | vVel | theta | salt | tile | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Float32 | Int64 | |
| 1 | 200 | 3600 | 187630.0 | 21119.0 | -1406.25 | 38.026 | 4.72379 | 3.0 | 0.114312 | -0.0385963 | 0.0101473 | 0.298643 | 35.0 | 1 |
| 2 | 200 | 7200 | 1.87491e5 | 21155.8 | -1406.25 | 37.9982 | 4.73115 | 3.0 | 0.114318 | -0.0387194 | 0.0102862 | 0.298643 | 35.0 | 1 |
| 3 | 200 | 10800 | 1.87351e5 | 21193.1 | -1406.25 | 37.9702 | 4.73861 | 3.0 | 0.114313 | -0.0388417 | 0.010443 | 0.298643 | 35.0 | 1 |
| 4 | 200 | 14400 | 1.87211e5 | 21230.9 | -1406.25 | 37.9422 | 4.74619 | 3.0 | 0.11431 | -0.0389631 | 0.0105994 | 0.298643 | 35.0 | 1 |
Super-impose trajectory over velocity field (first for u ...)
x=Γ["XG"].f[1][:,1]
y=Γ["YC"].f[1][1,:]
z=transpose(Γ["mskW"][1].*𝑃.u0);plt=contourf(x,y,z,c=:delta) plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
Super-impose trajectory over velocity field (... then for v)
x=Γ["XC"].f[1][:,1]
y=Γ["YG"].f[1][1,:]
z=transpose(Γ["mskW"][1].*𝑃.v0);plt=contourf(x,y,z,c=:delta) plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
6. Interpolate Velocities
dx=Γ["dx"]
uInit=[tmp[1,:lon];tmp[1,:lat]]./dx
nSteps=Int32(tmp[end,:time]/3600)-2
du=fill(0.0,2);Visualize and compare with actual grid point values – jumps on the tangential component are expected with linear scheme:
tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpx=fill(0.0,100)
for i=1:100
tmpx[i]=500.0 *i./dx
dxdt!(du,[tmpx[i];0.499./dx],𝑃,0.0)
tmpu[i]=du[1]
tmpv[i]=du[2]
end
plt=plot(tmpx,tmpu,label="u (interp)")
plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)")
plot!(tmpx,tmpv,label="v (interp)")
plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)")
And similarly in the other direction
tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpy=fill(0.0,100)
for i=1:100
tmpy[i]=500.0 *i./dx
dxdt!(du,[0.499./dx;tmpy[i]],𝑃,0.0)
tmpu[i]=du[1]
tmpv[i]=du[2]
endplt=plot(tmpx,tmpu,label="u (interp)") plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)") plot!(tmpx,tmpv,label="v (interp)") plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)")
Compare recomputed velocities with those from pkg/flt
nSteps=2998
tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);
tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);
refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);
for i=1:nSteps
dxy_dt_replay(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])
refu[i]=du[1]./dx
refv[i]=du[2]./dx
dxdt!(du,[tmp[i,:lon],tmp[i,:lat]]./dx,𝑃,Float64(tmp[i,:time]))
tmpu[i]=du[1]
tmpv[i]=du[2]
endplt=plot(tmpu,label="u") plot!(tmpv,label="v") plot!(refu,label="u (ref)") plot!(refv,label="v (ref)")
6. Compute Trajectories
Solve through time using OrdinaryDiffEq.jl with
dxdt!is the function computingd(position)/dtuInitis the initial conditionu @ tspan[1]tspanis the time interval𝑃are parameters fordxdt!Tsit5is the time-stepping schemereltolandabstolare tolerance parameters
tspan = (0.0,nSteps*3600.0)
#prob = ODEProblem(dxy_dt_replay,uInit,tspan,tmp)
prob = ODEProblem(dxdt!,uInit,tspan,𝑃)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
sol[1:4]retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 4-element Array{Float64,1}:
0.0
0.2007818920948153
2.2086008130429686
22.2867900225245
u: 4-element Array{Array{Float64,1},1}:
[37.525990625, 4.22379453125]
[37.52598925375626, 4.223794946058778]
[37.52597554129228, 4.223799094139329]
[37.5258384139952, 4.223840574222589]Compare recomputed trajectories with originals from MITgcm/pkg/flt
ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
maxLon=80*5.e3
maxLat=42*5.e3
#show(size(ref))
for i=1:nSteps-1
ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./dx;
plt=plot(sol[1,:],sol[2,:],linewidth=5,title="Using Recomputed Velocities",
xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false
plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,label="MITgcm Solution")
This page was generated using Literate.jl.